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Abstract 

The steady-state homogeneous vapor-to-hquid nucleation and the succeeding liquid droplet 



J> \ growth process are studied for water system by means of the coarse-grained molecular dy- 



namics simulations with the mW-model suggested originally in [Molinero, V.; Moore, E. B. 
J. Phys. Chem. B 2009, 113, 4008-4016]. The investigation covers the temperature range 
273 < T/K < 363 and the system's pressure ~ 1 atm. The thermodynamic integration 
scheme and the extended mean first passage time method as a tool to find the nucleation and 
^ . cluster growth characteristics are applied. The surface tension is numerically estimated and 

I is compared with the experimental data for the considered temperature range. We extract the 

nucleation characteristics such as the steady-state nucleation rate, the critical cluster size, the 
nucleation banier, the Zeldovich factor; perform the comparison with the other simulation 
results and test the treatment of the simulation results within the classical nucleation theory. 
We found that the liquid droplet growth is unsteady and follows the power law. At that, the 
growth laws exhibit the features unified for all the considered temperatures. The geometry of 
the nucleated droplets is also studied. 



*To whom correspondence should be addressed 



Introduction 



Nucleation is a fundamental process, which characterizes the mechanisms of the emergence of a 
new phase, and this is one of the most widespread ways, by which the phase transitions are initiated. 
Although a variety of theoretical descriptions for the nucleation exists, all of them are based on 
the same key idea: new phase starts to evolve within a mother phase from the nuclei, when they 
achieve such sizes and shapes, which facilitate the further stable growth of these nuclei. According 
to the classical nucleation theory (CNT), the stability of the nuclei is resulted by the confrontation 
of surface and bulk contributions in a free energy. This is relevant for the homogeneous scenario, 
which implies the equal probability for the appearance of a nucleation event over the sample, as 
well as for the heterogeneous scenario, where the some places in a sample are more attractive for 
the nucleation events (due to impurities, walls, etc.). 

Concerning the specific case of the homogeneous droplet nucleation during the vapor-to-liquid 
transition in water, there is the comprehensive experimental material due to series of investiga- 
tions (see, for example^"- and references therein). Here, the direct comparison of the experimental 
results with the predictions of the nucleation theories as well as with the data of the numerical sim- 
ulations performed by means of molecular dynamics (MD)^*^ and Monte Carlo— methods has 
revealed the noticeable discrepancies. So, for example, for the vapor-to-liquid nucleation in water 
at the identical conditions (pressure/supersaturation, temperature) the experiments, the theoretical 
models (CNT and others) and the numerical simulations yield the values of the steady-state nucle- 
ation rate Js, which differ by orders of magnitude. Under these circumstances, it could be quite 
reasonable to consider the features of the nucleation-growth kinetics in water at the molecular level, 
treating the vapor- to-liquid transition in the context of molecular interactions and movements. 

Recently, Molinero and Moore have suggested a coarse-grained "monatomic" model of water 
(mW), in which the anisotropy in the molecular interactions is simply realizing by means of an 
angular-dependent contribution.— The removal of the atomic interactions from the consideration 
accelerates the computations and, thereby, it inspires to probe the microscopic properties of the 
system on the extended time scales. Here, the phase transitions are convenient candidates to be 
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taken in handling. So, the homogeneous nucleation of ice was studied within the mW-model of 
water in Refs.— Therefore, it is tempting to extend these studies and to consider the details of 
the vapor-to-liquid phase transition on the basis of the mW-model. An important point is that the 
mW- model reproduces correctly the equation of state for the temperature range 250 <T /K < 350 
at /> ~ 1 atm. (see Fig. 4 in Ref.— ), which is relevant at the consideration of the droplet nucleation 
in water. 

From viewpoint of the CNT, three principal parameters are enough to restore the basic aspects 
of the steady-state nucleation. These can be, for example, the steady-state nucleation rate Js, 
the nucleation barrier AG and the Zeldovich factor Z. Of course, those three parameters can be 
taken in another combination (for example, the "reduced moment", the lag-time, and the steady- 
state nucleation rate, like it is suggested in Ref.—). Nevertheless, the surface tension o, which 
characterizes the interphase layer and contributes to the nucleation barrier through the surface 
free energy term, requires the independent treatment. In the direct computer simulations, the 
different adapted convenient approaches based on the Fowler formula, the Kirkwood-Buff formula 
and others are utilizing to define accurately the surface tension.-i^*22 However, there is a necessity 
at the study of nucleation to apply such a method (i) that gives a possibility to estimate the surface 
tension from the raw simulation data, (ii) that is applicable to characterize the surfaces of the 
microscopical nuclei with a pronounced inherent curvature, and (iii) that considers the genuine 
interphase (vapor-liquid) properties without reference to a vacuum phase. 

In the present work, we study the nucleation-growth processes of water droplets on the basis 
of MD simulations with the mW-model. To define the parameters of the nucleation and the droplet 
growth, we apply the statistical treatment of the simulation data on the basis of the thermodynamic 
integration scheme and the mean first passage time (MFPT) approach. Similar to the thermo- 
dynamic integration scheme, the MFPT approach utilizes the time-dependent configurations as 
resulted from the independent runs under identical conditions, however, the MFPT is focused on 
the averaged time scales, at which a system characteristic (reaction coordinate, order parameter) 
appears for the first time.—"— We show that the thermodynamic integration scheme and the MFPT 
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method provide a convenient tool to treat the simulation results (and/or the experimental data) con- 
cerning both the nucleation and the growth kinetics. For the considered case of water, we define the 
set of the characteristics for steady-state homogeneous nucleation and growth of the liquid droplets 
on the basis of MD simulation data. 

Numerical schemes 

Thermodynamic integration. - The surface energy w can be defined as an excess energy per unit 
area of the surface that is conditioned by the lack of neighbors for the surface particles in compar- 
ison with the bulk particles (see Ref.— ). If one restricts the consideration by the closest neighbors 
only with the pairwise additive interactions u{rij), then the following relation appears directly 

w=^u{rij){z-z)n, (1) 

where r'jj is the average distance between the neighbors in a new phase, the quantity n' denotes the 
number of surface particles per unit area and depends on the size of a nucleus, z and z' are the first 
coordination number of bulk and surface particles, respectively. Then, the surface tension can be 
estimated directly by the thermodynamic integration of the surface energy as 

The reaction coordinate A or the so-called A -scaling— is associated with the rescaled cluster size, 
X = {n/n*y/^ , which is equal to zero if there are no nuclei in the system and to unity if the nucleus 
size has the critical value n*. The notation {■ ■■)x means an ensemble average at a particular value 
of A. 

Extended mean first passage time method. - According to the continuous Zeldovich-Frenkel 
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scheme, the nucleation process can be described within a Fokker-Planck-type equation 



dt 



Nnjt) 



(3) 



where n is the cluster size, Nn{t) is the time-dependent cluster size distribution over unit volume, 
Jn is the current over cluster size space, g+ is the monomer attachment rate to a n-sized cluster and 
Nn"^ = Nq' exp{—l5AGn) is the equilibrium cluster size distribution, AG,, is the work required to 
form the n-sized cluster and /3 = 1/ {ksT). 

If one considers the n-dependent term AG„, the nucleation regime is directly associated with 
the vicinity of critical value of the cluster size, n* , where the term AG,,* corresponds to a nucleation 
barrier and has a maximum. Assuming that the nucleation barrier can be expanded into the Taylor 
series in this vicinity 



AG,, = AG,,* + £ 

k=2 



the approximated evaluation of Eq. (??) in the vicinity of nucleation regime can be written as 



exp(/3AG,,* 



X / dn exp 
/o 



(5) 



k^2 



kl 



dn'' 



The series in the exponential of Eq. Q contains an information about the geometrical peculiarities 
of the term AG,, around its maximum at n*. Namely, the second contribution of the series is related 
with the Zeldovich factor Z and characterizes the curvature of the barrier at the top 



2 dn^ 



(6) 



d) (2) 

Moreover, the ratio of the third and the second contributions, which is AG)jJ^*/3AG)j^,j*, indicates 
on the asymmetric properties of the barrier. For example, if the ratio is equal to zero, then the 
barrier is symmetric one and can be approximated by a parabolic geometry. This means for the 
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given example that we are restricted here only by a case with k = 2, which corresponds to the 
Zeldovich approximation. Here, the analytical expression for the steady-state nucleation rate Js 
can be directly obtained from Eq. ([5]) within the MFPT method,— where the averaged time scale 
of the first appearance of the n-sized cluster x^^^^ considered: 

^''^ = ^{l+erf[v/^Z(n-n*)]} (7) 
^ erfc[v^Z(n-n*)]. 



2JsV 



Here, V is the system volume, and erf(;c) = 2k^^/^ [q exp{—t^)dt is the error function. 

The MFPT method provides the next useful capabilities in the treatment of the nucleation- 
growth processes. The first one is related with the critical value n*, which is located at the inflection 
point, i.e. at the point, where the first derivative ((9t^^^^/ dn)n=n* has a maximum. Thus, a simple 
analysis of dx^^^^ j dn yields the critical value n* (see Fig. [I]). For the particular case of Eq. (|7]), 
one obtains directly that n = n*, when T,^,f ^ = 1/ (2^^^) that is the consequence of the nucleation 
barrier symmetry. The second property is that the Zeldovich factor can be directly extracted from 
MFPT as 

Z=J.V^ . (8) 



dn 

The geometric constructions, corresponding to this equation, are presented in Fig.[TJ Equation (??) 
indicates that the smaller values of Z are resulted from the smaller values of {dxj^^^^ / dn)n=n* at 
the fixed JsV . On the other hand, the smaller values of the Zeldovich factor correspond to the flatter 
nucleation barrier curve AG,, near the critical size n* . And, finally, the third property is associated 
with the steady-state nucleation rate Js, which can be defined from the MFPT distribution as Js = 
I / (x^f^'^^V) at n, where {dx^^^^ / dn)n>n* approaches the minimum and the distribution x^^^^ 
starts itself to demonstrate a steady-like n-dependence (see Fig. [T]). Thus, using the known mean 
first passage time distribution x^^^^ one can directly define the critical value n*, the Zeldovich 
factor Z and the steady-state nucleation rate Js by a direct numerical analysis. 

Nucleation-growth kinetics is characterized by the nucleation time scale x„ = 1/ {JsV) and the 




Figure 1: (Color online) Top: Schematic plot of the MFPT distribution for the cluster size n as ob- 
tained from simulation (or experimental) data. The regions I and II are associated with nucleation 
and cluster-growth regimes, respectively. The routine for finding the nucleation characteristics 
from the MFPT-curve is presented. The gentle slope of the MFPT-curve at the transition value, 
n = n*,is evidence of the smooth form of the nucleation barrier AG„ in the vicinity of n*, while the 
location of the inflection point (full circle) below the half-height 1 / {2JsV) indicates qualitatively 
on the barrier asymmetry. The pronounced increase of MFPT-curve in the region II appears due to 
the fact that nucleation and cluster-growth time scales are comparable, and this part of the curve 
as an inverted one, n(T^^^^), can be used to estimate the parameters of cluster growth kinetics. 
Bottom: First derivative of the MFPT distribution, dx^^^^ /dn. Here, the maximum is associated 
with the inflection point, which is directly located at the critical value of cluster size, n = n*. Po- 
sition of the next extremum (minimum) on {dx^^^^ / dn)n>n* corresponds to the nucleation time 
scale Xn = l/{2JsV) as defined from the main MFPT distribution x^^^^ . Inset: Typical cluster 
growth curves obtained from the independent simulations. 

cluster-growth time scale Xgr. The ratio between these time scales distinguishes the separate cases 
for the numerical treatment within the MFPT method: (i) If T„ ^ Xgr, then x^^^^ demonstrates a 
clear defined plateau of the height T„ = l/{JsV), that simplifies significantly accurate estimation 
of the nucleation rate; (ii) If these time scales are comparable, T„ ~ Xgr, then the errors can appear 
in the estimation of x„, since the boundary between nucleation and growth in MFPT distribution is 
smeared. 

Furthermore, the MFPT method gives a convenient tool to extract the characteristics of nucleus 
growth kinetics, which follows the nucleation regime in the MFPT distribution (see Fig.[T]). In fact. 
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the inverted MFPT distribution, «(t^^^^), has the statistical meaning of the most probable cluster 
growth law for the growth regime of the MFPT curve. 

Following Ref.,— the growth law of a cluster can be taken in general form as 

R{t)=R, + {%it)\ (9) 

where R and R^ is the radius of the growing cluster and the critically-sized cluster, respectively; v is 
the growth exponent and is the growth constant, which has a dimension of [/m^/^/^]. Then, the 
growth rate is G{t) = v'i^^^t'^^^ , while the acceleration of a cluster growth can be formally defined 
as a{t) = v(v — l)'^^^it^^^. The steady cluster growth with a constant growth rate corresponds to 
the particular case of v = 1, where the growth rate coincides with the growth factor, i.e. G{t) = 
= const, otherwise (at v 7^ 1) one has the process with unsteady growth rate.— Further, taking 
into account that the volume of a growing cluster evolves with time as V{t) = Cg[R{t)]^ and N{t) = 
pc.iV{t), where Cg is a dimensionless cluster-shape factor (cg = 4;r/3 in a case of the sphere) and 
pel is the density of the cluster-phase, one can write the growth law in the extended form: 



n[t,tcj 



1 + it - r^^P^ + 3^,f (r - r.r (^) ' + 3^,^^ " ^c) ^ (^) ' 

'(10) 

Here, the lag-time tc defines the appearance of the critically-sized cluster. Then, the term n{T^^^^) 
can be fitted for the growth regime by Eq. (??) to extract the growth characteristics: the cluster- 
shape factor Cg, the growth constant '^ci and the growth exponent v. At rapid growth of small 
clusters the last two contributions in Eq. (??) can be neglected, and the growth law takes the 
form— 

n{t,t,)c^n*+cgp,0^^{t-t,)^\ (11) 

where V is positive. 
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Figure 2: (Color online) Top: Direct MFPT distributions (stepwise-curves) and their interpolations 
(smooth solid curves) for the temperatures T = 273, 293 and 333 K. Dots on the curves indicate the 
inflection points, which define the critical sizes n*, and the time scales T = (/jV)"^ corresponding 
to the nucleation rates. Thick short lines are the linear parts of the interpolated curves near n* and 
define the ranges of errors in critical sizes. Note that the errors in the nucleation rates can be also 
defined as a result of changes in t = (JsV)^^ due to the correction of different interpolations with 
the same accuracy in the reproduction of the direct MFPT-distributions. Bottom: First derivative 
of the MFPT distributions, dx^^^^ / dn, for the same temperatures. 

Computational details 

Molecular dynamics simulations were performed in the spirit of previous studies of the structural 
transformations in this system described in Refs.— »^ with only difference in the details related with 
the considered thermodynamic range. We have examined the system composed N = S 000 particles 
(molecules) interacting via the mW-potential in the cubic cell with the periodic boundary condi- 
tions in all directions. The time-step for numerical integration was 1 fs; and the NpT (number, 
pressure, temperature) ensemble was applied with p = I atm. Pressure and temperature were con- 
trolled via the Nose-Hoover barostat and thermostat, respectively, acting uniformly throughout the 
system. The damping thermostat and barostat constants were taken to be Tj- = Tp = 10 fs. The 
parameters of the mW-potential are completely identical to those reported in Refs.— »^ 
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Initially, the set of a hundred of independent samples was prepared and equilibrated at the 
temperature T — 900 K on the time scale 50 ps (i.e. 50 000 time-steps). The correspondence of 
the systems to the vapor phase was directly confirmed by the particle diffusivity and the distinc- 
tive particle radial distribution functions. Moreover, following Ref.,— the samples were cooled at 
10 K/ns to the desired temperatures from the range 273 <T /K < 373 (at /> ~ 1 atm.)y Then, over 
a time scale ~ 1 ^ 10 ps each a system was 'equilibrated' till the disappearance of the pronounced 
fluctuations in temperature and pressure, after that the initial configurations were stored for the 
further study of the vapor-to-liquid nucleation process. Note that this cooling procedure is similar 
to the reported one in Ref. The following A^pT-simulations starting from these configurations - 
a hundred for each considered temperature - were performed to collect the statistics of the inde- 
pendent nucleated events, where the time-dependent cluster size distributions Nn{t) were evaluated 
(for an every run). The averaged time scale for the simulations in this nucleation-growth regime 
was 50 ns. On the basis of the found A',j(?)-distributions, the MFPT-curves were extracted and the 
nucleation characteristics were estimated according to the scheme presented above. After this, the 
critical sizes n* defined from the MFPT-curves were used at the retreatment of the simulation data 
with the aim to define the distributions of the energy (O over the reaction coordinate A . 

An identification of the particles, which belong to liquid phase, was performed in the spirit of 
the Stillinger rule.^^ First, the particles are "neighbors" (or bonded) if the distance between their 
centers is less than r^, where is the position of the first minimum in the pair correlation function 
of the liquid phase (at the same conditions). Further, a particle is considered as a liquid-like if it 
has, at least, four neighbors 0. 

'it is necessary to note that the mW-model reproduces correctly equation of state p{T) for this temperature range 
(see Fig. 4 in Ref.ii). 

^The last condition allows one to remove from the consideration those particle-pairs, which are result of the instant 
random event and are not related to the formation of a new phase. 
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Figure 3: (Color online) (a) Radial distribution function g{r) of liquid water as resulted from the 
mW- model: results for bulk water at T = 298 K and p = (full circles) reported in Ref., data 
for the bulk range of critically- sized droplets at T = 293 K and p = I atm. Inset: Distribution 
of the first coordination number for the water molecules of a critically-sized liquid droplet at the 
temperature T = 293 K. The histogram corresponds to the total distribution; the line with triangles 
present an impact from the bulk molecules, z; and the line with rotated triangles shows the contri- 
bution of the surface molecules, z'. The data are averaged over set of runs, (b) Main: Temperature 
dependence of the surface tension o. The simulation results show the averages (full circles) and 
standard deviations (error bars) from independent runs; experimental data are presented by open 
circles, whereas the dotted line is the interpolation by o{T) =B[{Tc - T) jT^^i 1 ^h\{Jc-T) /T'c] }, 
B = 235.6 N/m, b = -0.625, m = 1.256 and Tc = 647.15 K.-^ Inset: Surface energy and slope of 
the surface energy d(o/dX as functions of A = {N /NcY^'^ for a growing water droplet. At the 
critical size N — Nc one has A = 1. The presented results are outcome of a single simulation run. 

Results 

In Fig. |2] the obtained MFPT distributions x^^^^ and their derivatives dx^^^^ /dn at the temper- 
atures T = 273, 293 and 333 K are presented as an example. On the basis of the defined values 
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of the critical size n* and nucleation rate Js = {t V)^ , the values of the Zeldovich factor 
were extracted within Eq. (??). The typical A -dependences of the surface energy and its deriva- 
tive as obtained from a single run are shown in the inset of Fig. Sb). The averages of the slope 
{d(o/dX)i over different runs were used to estimate Eq. (??) by means of the trapezoidal method. 
The smooth character of the curves allows one to restrict oneself by the method and to exclude 
higher order integration schemes.— It is necessary to note that errors in the critical size n* have 
not been considered at the estimation of the surface tension with Eq. (??). 

Coordination number. - The inset of Fig. (S^a) shows the distribution of the first coordination 
number for the water molecules generated a droplet of the critical size in the system at the temper- 
ature T = 293 K. The presented histogram is the cumulative result of the distributions for surface 
and bulk molecules. Remarkably, these distributions (for bulk and surface molecules) are sym- 
metric ones as well as reproducible by the Gaussian functions. Further, the averaged values of the 
coordination numbers z and z' are extracted from the distributions and appear to be 5.8 and 3.93, 
respectively, for the case. Other important observation is that the term z is practically unchange- 
able with temperature, whereas the coordination number in a surface layer demonstrates a smooth 
insignificant decrease with the temperature increasing [for comparison, from z'{T = 293 K) = 3.93 
toz'(r = 353 K) = 3.47]. 

Moreover, the value of the coordination number for bulk molecules in the droplets coincides 
with the found value of z extracted on the basis of the integral definition— 



where rc is the first minimum position in the radial distribution function g{r). Nevertheless, the 
found value z{T = 293 K) = 5.8 differ from the result of Molinero and Moore^ obtained within 
the mW-model for the bulk water, z{T = 298 K) =5.1. To understand the reasons of the discrep- 
ancy, we compare the corresponding radial distribution functions in Fig. HJ^a). As can be seen, 
the intensity of the first maximum of g{r) is higher for the case of the water droplets, although the 




(12) 
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maximum is located at a lower distance. This feature indicates that the liquid phase is characterized 
by the more pronounced short-range ordering for the case of the microscopically small nucleated 
clusters than for the equilibrium liquid phase considered in Ref.— 

Surface tension. - In Fig. Ob), the temperature dependence of the surface tension of the 
critically-sized nuclei is presented. The results obtained from simulation data demonstrate the 
known decrease of the surface tension with temperature □ and reproduce precisely the experimen- 
tal data of this term for a planar liquid-vapor interface.— Such a good agreement of our results 
with experimental data is unexpected because of the two next reasons, mainly. First, in contrast 
to an inherent water system, the mW-model excludes the long-range intermolecular interactions, 
which still can have an influence on the interface effects. Second, no adjustment of simulation data 
was performed to take into account the finite size effects. Thereby, the corrections to surface free 
energy in the spirit of the Tolman's ansatz, which are appeared to be proportional to the inverse 
linear size of the critical nucleus, were unconsidered for the surface tension. 

On the other hand, the surface tension of a planar interface was recently defined with the mW- 
model for the temperatures 250 <T/K< 350.— The values of the surface tension reported in 
Ref.— have the lower values in comparison with the values presented in Fig. [2b) for the water 
droplets, and the difference is about 5^10 percents for the same temperature range. Although this 
difference could be attributed to the Tolman length with the negative value, like it was reported by 
Kiselev and Ely— for the liquid-ice surface tension, we assume that our values of o can be over- 
estimated because of the neglect the three-particle interactions by computational protocol within 
Eq. (??). 

Nevertheless, the surface tension decrease with the temperature is directly consistent with the 
temperature decreasing of uijij), of the surface coordination number [or the increase of the differ- 
ence (z — z')] and of the surface particle density, where the last two contributions are practically 
counterbalanced by each other^. 

^According to results of Ref. '''^^ obtained for the planar surface tension at a single temperature T = 300 K, the 
mW-model gives the best agreement with the experiment in comparison to the models: SPC, SPC/E, TIP3P, TIP4P, 
TIP5P 

^Assuming the spherical droplet of the radius R with the thickness of the surface layer A, the surface particle density 
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Recalling the previous debates,— the temperature range investigated here can contain the in- 
flection points in the vapor-liquid surface tension of water. The results, presented in Fig. Hb), 
indicate on the absence of the clear detected inflections in the T-dependence of the surface ten- 



sion. 
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Figure 4: (Color online) Temperature dependence of the homogeneous droplet (vapor-to-liquid) 
nucleation rate Js{p — l atm, T) in water. Comparison of the simulation results with the mW-model 
[the vapor density is Pv G [1.12; 1.55] x 10~^nm^^], with the atomistic SPC/E-model from Ref.— 
[the density is p,, G [1.23; 1.86] x lO^^nm"^], with the atomistic TIP4P-model from Ref.^^ [the 
numerical density is p,, = 1.55 x lO^^nm^^] and the treatment of nucleation data within the CNT. 
Solid line show the best fit by means of the function Js{p — 1 atm, T) oc exp(— 0.00092 T^-^). 

Steady-state nucleation rate. - In Fig. H the steady- state nucleation rates obtained from the 
MFPT treatment of the simulation data with the mW-model are presented as a function of tem- 
perature (see also Table[T]). The presented data cover the density range p,, from 1.14 x lO^^nm^^ 
to 1.55 X lO^^nm^-^. As can be seen, the values of Js{T) compare well with those obtained for 
the same density range by Matsubara et al. with the atomistic SPC/E-model.— Nevertheless, as 
contrasted to results of Ref.,— which are scattered over {Js, r)-plot, the values of the nucleation 
rate Js obtained within the mW-model demonstrate the smooth decrease with temperature over 
can be roughly estimated n'{R) ^[R — {R—A)^/R^]. 
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the considered temperature range. Moreover, this decrease is well-reproduced by the dependence 
ln[7,(r)] = -0.00092 T^-^ + 75. Among all the data presented on Fig. gfa), the highest value of 
the nucleation rate appears from the simulations of Ref.— with the TIP4P-model. 

On the other hand, it is attractive to test the dependence Js{T) within the CNT treatment with 
the extracted values of the other nucleation characteristics. So, the original Becker-Doring formu- 
lation yields 

7f^^ = ^./^expf-^y (13) 



pi V Ttm \ kgT 
where the barrier can be taken as 

^ = 37t{n*Z)^ (14) 

and (Too is the surface tension for a planar liquid-vapor interface.^^ As can be seen from Fig. |4l 
although the both dependencies demonstrate a similar behavior decaying with temperature, the 
pure simulation results for nucleation rate have in two orders higher values in comparison with 
J^^^ . This is evidence of the difficulties at the description of the droplet nucleation in water by 
means of the CNT, which are similar with those reported earlier for the studies with the atomistic 
models^^*^ as well as with the experimental data (see Fig. 3 of Ref.—). 

Critical cluster (droplet) size. - Figure |5] illustrates the temperature dependence of the critical 
cluster size n*, where the simulation results with the mW-model are compared with the simulation 
data of Matsubara et ah— and of Yasuoka et al.^ as well as with the predictions of the Kelvin 
equation 

n — ; — ; — rrr- (15) 

3 pf[kBT\n{p/p^)]^ 

Here p^ is the saturated water vapor pressure.— 

First, for the mW-model the critical cluster size reveals a slight decrease with temperature 
from n* = 75 to 40 particles over the temperature range 273 < T /K < 363. This change of the 
cluster size means the decrease of the droplet radius from 4 to 3.3 of the averaged water molecule 
diameters. Obviously, the change is insignificant. Moreover, the observed decrease is masked by 
errors, which were defined as the curvature range width in the MFPT-distributions (see Fig.|2l). The 
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Figure 5: (Color online) Temperature dependence of the critical cluster size, determined from sim- 
ulations within the different models of potential fields (mW-model, SPC/E, TIP4P— ), compared 
to prediction by Kelvin equation. In the case of the mW-model, error bars are defined by a width 
of the curvature range in MFPT-distributions. 

range of errors is ±(10^25) particles, that is awaited to be reasonable, since it covers only a few 
surface part of a water droplet (Fig. [6l). In addition, according to the definition within the MFPT- 
method these errors should be considered as the probable deviations from n* in a statistical sense. 
The comparison with the results obtained for the atomistic models (TIP4P and SPC/E) reveals that 
the values of n* obtained within the mW-model overestimate the data of the SPC/E-model, but are 
in agreement with a single value found by Yasuoka et al. in simulations with the TIP4P-model. 

Further, the mW-model result for n*(r)-curve is different from the predictions of Eq. (??), 
which yields the increase of n* with the temperature T (see Fig.[5])y. However, it should be noted 
that as far as the predictions of the Kelvin equation are concerned, it gives the values n ^ 2 ^ 12 
molecules for the temperature range 273 <T/K< 363 and the pressure p = I atm. In terms of 
the linear cluster sizes, these values correspond to 1 ^ 2 water molecule diameters. It is clear that 

^The evaluation of the supersaturation S = pyj p\ was performed with the experimental values of the saturated 
water vapor pressure p\,. Nevertheless, one needs to note that the mW-model can yields the results different from the 
experimental data for p'l{T). 
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14.1 nm 



2.07 nm 



Figure 6: (Color online) Left: Snapshot of the water system at the temperature T = 293 K and 
the pressure p = I atm. at the moment, when the critically-sized droplets are appearing. Right: 
Growing droplet of the same system. 

the treatment of the stability of such a small cluster from the thermodynamic point of view, which 
requires the availability of the separated surface and bulk regions of the cluster, is impossible. ^ 

Nucleation barrier and the Zeldovich factor - The temperature dependence of the next nucle- 
ation characteristic, the Zeldovich factor Z, is presented in Fig. |7J As can be seen, this quantity 
decreases from the value 0.028 to 0.014 with the decrease of the temperature T (see also Table[T])- 
Such a behavior is a direct evidence that the nucleation barrier loses its sharpness and becomes 
more smoother with the decreasing the temperature. This is qualitatively in an agreement with 
the prediction of the CNT.— Moreover, if one assumes that the CNT yields the correct results for 
the vapor-to-liquid nucleation of water then it is possible to define the nucleation barrier by means 
of the simple relation (??). The direct evaluation yields the correct tendency of the temperature 
dependence for the nucleation barrier, though this tendency is not so pronounced one as it could 
be expected for a sufficiently wide temperature range considered here. According to the CNT, 
the nucleation driving force |A/i|, growing with the decrease of T, must reduce the nucleation 
barrier.— 

Within Eq. (??) and simulation results one has that the nucleation barrier decreases from 



^ We note that the direct comparison of the predictions of the Kelvin equation with simulation results should be 
considered as very approximate, since the saturation curve resulted from a considered model can be different from the 
real water saturation curve.— For the mW-model, the additional studies are necessary to clarify this point. 
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pAG^*^ = 12 ± 2.7 to pAG^.^ = 9.7 ± 2.5 with the decrease of the temperature from T = 363 
to 273 K. On the other hand, the temperature dependence of the nucleation barrier as predicted by 
the CNT is defined by 

AG„* 



pf\A^r 

where |A/i| is the chemical potential difference of particles in the vapor and in the liquid phase. 
So, the observed behavior of the nucleation barrier can be explained for the case, where the change 
of (J^/pf with the temperature is completely counterbalanced by the change of |A/ip. 

Remarkably, the comparable values for the nucleation barrier arise with the atomistic SPC/E 
model (see Table I in Ref.— ), where the barrier changes from /3AG„* =8.1 to /3AG„* = 6.7 with 
the temperature decreasing from T = 325 to 275 K. Thus, the observed results for the nucleation 
barrier can not be considered as a consequence of the coarse-grained character of particle interac- 
tions in the mW-model. 

Table 1: Simulation results: system temperature T (K); vapor number density p,, (x 10^^ nm^^); 
critical cluster size n*; nucleation barrier AG/hgT; nucleation rate Js (x 10-^^ m^-^s^^); the Zel- 
dovich factor Z. 



T 


Pv 


nc 


AG/ksT 


Js 


Z 


273 


1.548 ±0.095 


75 ±25 


9.72 ±2.55 


0.35 


0.0135 ±0.0009 


283 


1.462 ±0.060 


71 ±27 


9.76 ±2.06 


0.30 


0.0143 ±0.0008 


293 


1.426 ±0.077 


65 ±27 


9.78±2.39 


0.26 


0.0157 ±0.0008 


303 


1.425 ±0.051 


58 ±22 


10.07 ±1.57 


0.24 


0.0182 ±0.0007 


313 


1.353 ±0.062 


55±21 


10.21 ±2.25 


0.21 


0.0189 ±0.0009 


323 


1.313 ±0.047 


52±18 


10.43 ±2.18 


0.19 


0.0202 ±0.0008 


333 


1.250 ±0.037 


50±15 


10.98 ±2.53 


0.17 


0.0216 ±0.001 


343 


1.224 ±0.032 


45 ±13 


11.10±1.70 


0.15 


0.0241 ±0.0007 


353 


1.174 ±0.045 


42 ±12 


11.50±1.94 


0.13 


0.0262 ±0.0009 


363 


1.140 ±0.033 


41 ±12 


11.97±2.17 


0.11 


0.0275 ±0.001 



Growth laws of the nucleated droplets. - Bottom inset of Fig. |7] shows the growth curves of 
the liquid droplet at the different temperatures, which were found from the statistical treatment of 
simulation data by means of the MFPT approach as it was discussed above. Hence, these curves 
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Figure 7: (Color online) Main: Temperature dependence of the Zeldovich factor Z as defined by 
the MFPT method on the basis of the simulations with the mW-model. Top inset: Temperature 
dependence of the growth factor which is found from the fit of Eq. (??) to the simulation data. 
The growth exponent v = 1.3 and the term A = c^pii^Scf"^ In' = 1.16 ±0.2 appear to be invariant 
respective the temperature. Error bars show the standard deviations from the averages. Bottom 
inset: Growth curves of the liquid water droplets emerging in the vapor phase at the temperatures 
273 <r/i^< 373. 

depict the most probable growth laws in a statistical sense. As can be seen from the figure, at 
the lower temperature the droplet growth occurs faster. At the same time, all the curves are well 
reproduced by Eq. (??), and the fitting to the simulation data yields the following features. The 
growth exponent v in Eq. (??) appears to be invariant over the temperature T and it takes the value 
V = 1.3 at all the considered temperatures. This indicates that the increasing the linear size, i.e. 
the radius, averaged over all the directions of the liquid droplet follows for all the temperatures the 
growth law R(t) = (?fc/0''^' herewith, the droplet growth itself is unsteady. 

Further, the growth factor ^f^/ decreases with the increase of the temperature T (see top inset 
of Fig.|7]), that characterizes the faster droplet growth at the lower temperatures. Following Ref.,— 
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Eq. (??) can be written in the rescaled form: 



n(<^)/n*~l+A(^-l) 



3v 



(16) 



where A = CgPi{^cltc) /n* and ^ = t/tc is the rescaled time. We found that in accordance with 
the rescaled form (??) all the growth curves collapse onto a single curve independently of the 
temperature. In addition, the parameter A in Eq. (??) takes the same value for all the considered 
temperatures, i.e. A = 1 . 16 ± 0.2. This can be evidence of the generic features of the water droplet 
growth process. Remarkably, this result is correlated with the features of the crystal growth kinetics 
for a model glassy system under shear drive, which were reported in Ref.— 

Shape and sphericity of the nucleated droplets. - Other issue, which is crucial in the CNT, 
is related with the shape and the anisotropy of the growing droplets .-^^^ A convenient way to 
perform this study in our case is to use the asphericity parameter in the next definition: 



is the components of the moment of inertia tensor associated with a droplet, mo is the molecule 
mass, CU,j8 e {x,y,z} are the components of the vector r between the droplet center-of-mass and 
molecule z; the brackets (...) mean the statistical average over critically-sized droplets of the dif- 
ferent simulation runs. This definition of a asphericity parameter sets the characterization: for 
a spherical droplet one has Sq = 0, whereas for an elongated and string-like cluster one obtains 
5o — 1 . We found that independently on the particular conditions (temperature, vapor density) the 
asphericity parameter for the considered (p, r)-range is Sq ~ 0.008 ±0.0002. It indicates on the 
nucleated water droplets of the sphere-like form, which is also confirmed by a visual inspection 
of snapshots (Fig. [6]). We remark here, this result is not the same with the findings of Refs.,— »^ 




where 



lap =Y.'^o{rf5ap - 



7=1 
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where the detected critical clusters in water had the significant deviations from a spherical form. 
A possible reason affecting the observed discrepancy could be different cluster definitions applied 
by Matsubara et al. in Ref.— and used in the present study within the statistical treatment. In 
addition, the low values for the nucleated droplets were obtained in Ref.,— ~ 16 ^ 22 partilces, 
and the pronounced deviation from a spherical form can be considered as a signature of the finite 
size effects: a weak structural rearrangement in such a system raises the significant change of its 
shape. 

Conclusions 

The coarse-grained models for particle interactions in molecular systems provide the good oppor- 
tunity to study early stages of the phase transitions by means of the numerical simulations. In this 
work, the processes of the steady-state homogeneous vapor-to-liquid nucleation and the growth 
of liquid droplets in water were considered within the mW-model, which treats the molecular 
interactions excluding any details of the direct oxygen-hydrogen interactions and electrostatics. 
Despite the apparent coarsening in the description of the molecular interactions, we have shown 
that the mW-model provides interesting information concerning the droplet nucleation in water va- 
por, thereby complementing the simulation results obtained earlier within all- atom models of water 
such as TIP4P and SPC/E.— '^>^>^ It is necessary to note, the results reported here are obtained on 
the basis of the extended statistical treatment within the MFPT approach and the thermodynamic 
integration scheme. 

The surface tension of the nucleated droplets was computed within an approximation that is 
restricted by the consideration of the two-particle interactions only without handling the three- 
particle contribution to the energy of system. The obtained values demonstrate the decrease of 
the surface tension with the temperature growth. It is necessary to note that the applied numerical 
scheme gives the higher values for the liquid- vapor surface tension of the droplets in comparison 
with the values for the liquid-vacuum surface tension of a planar interface reported in Ref .— We 
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suppose that this difference is rather a result of the approximations applied to the surface tension 
definition than the mW-model product. 

Further, the evaluated values of the steady-state nucleation rate are comparable with the results 
for all-atom models as well as with the treatments within the classical nucleation theory. Unfor- 
tunately, we could not to perform the direct comparison of the obtained nucleation rates with the 
experimental data, because we found no experimental Js for the {p, r)-line considered here. Nev- 
ertheless, quantitative extrapolation of the obtained outcomes indicates on the difference between 
simulated and experimental results, that is similar with the known difference between the experi- 
mental data and the CNT predictions. -^i^ On the other hand, for the critical size of the nucleated 
droplet we found the values within the range 30 ^ 100 particles, which are expected to be compa- 
rable with the experimental data (Fig. 13 in Ref.-). So, the additional studies are highly desirable 
in this field. 

According to our results, the growth of nucleated droplets in the system is characterized by the 
remarkable features: the growth law of the droplet radius follows the power law, R{t) t^^^, and 
the growth is not steady (with the time-dependent growth rate G{t) oc t^-^). Moreover, the simple 
rescaling on the critical droplet characteristics yields the unified form of the growth law at all the 
considered temperatures. 

Finally, we found that the critically- sized droplets have a shape, which is close to spherical one. 
Note, that the deviations from spherical shape of the water droplets at homogeneous nucleation, 
which were established for the SPC/E-model by Matsubara et al. (see Ref.—), could be simply 
originated from the extremely low obtained values for the critical size n*. 
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